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FOREWORD 


The  research  reported  herein  was  perfoimed  by  North 
American  Aviation,  Inc.,  Columbus,  Ohio  for  the  Hypersonic 
Flight  Section,  Flight  Branch  of  the  Flight  Dynamics  Labora¬ 
tory,  Wright  Air  Development  Division.  The  work  was  accom¬ 
plished  under  Air  Force  Contract  No.  AF  33(6l6)~7100, 

Project  No.  1366,  Task  No.  13967,  "A  Study  of  Aerodynamically 
Oriented  and  Stabilized  Satellites . "  This  research  was 
carried  out  by  the  Engineering  Research  and  Aerothermodynamics 
Development  Groups  of  the  Columbus  Division,  North  American 
Aviation,  Inc.,  with  Dr.  D.  M.  Schrello,  Engineering  Research 
Group,  as  Project  Engineer.  Mr.  Joseph  Ondrejka,  Flight 
Dynamics  Laboratory,  was  WADD  Project  Engineer. 

The  results  of  this  study  are  reported  in  a  series  of 
three  volumes  of  which  this  is  Volume  III.  The  other  reports 
in  this  series  are: 

VOLUME  I:  "Librations  due  to  Combined  Aerodynamic 

and  Gravitational  Torques,"  by  D.  M. 

Schrello 

VOLUME  II:  "Aerodynamic  Analysis,"  by  Paul  H.  Davison 
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ABSTRACT 


A  computer  program  is  presented  which  integrates  numerically 
the  pitch  equation  of  an  aero dynamically  stabilized  satellite. 

The  mathematical  theory  of  the  linear  second  order  equation  with 
periodic  coefficients  and  the  integration  procedure  are  dis¬ 
cussed  to  clarify  the  structure  and  application  of  the  program. 
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INTRODUCTION 


The  equation  of  satellite  pitching  motion  derived  in  Volume  I  of  this 
series  of  reports  is : 

0"  ♦  l  Cl  8’  ♦  {  3M  [  I  -  e  v!  *  Pi  9  -  2  a  +  YP  (l) 

where  V  ,  the  true  anomaly,  is  the  independent  variable, 

9  is  the  pitch  angle  from  the  local  horizontal, 

£  is  the  orbital  eccentricity, 

M  is  the  inertia  ratio  parameter, 

and  the  remaining  quantities  are  defined  by  the  following  relations: 

r*  .  (2) 

/♦  6  e*t.  v 


P-  0>r  (/>//>,)(  Vz/ \/)X  [  1+  26  (I-  C+4V-)]  , 

Q.  *  Up  (y©/ f>r  )  (  Vh  /  V)  [  I  -  £  ( I  -  out v )  ]  ~  £  a**)  xr  , 


in  which 

(Vfi/V)  *  /  -  (  u?£  /  i  , 

<Tp  =  ±  /  /u  (/»*)"  ,  <3) 

r,  /  /> 

the  subscript  P  denotes  perigee  conditions, 

JJL  is  the  product  of  the  universal  gravitational  constant 
and  the  mass  of  the  earth, 

U>e  is  the  Earth’s  rotation  rate, 

i  is  the  orbit  inclination, 

p  ,  the  atmospheric  density,  is  a  tabular  function  of  the 
/  altitude  h  ,  see  for  instance  Minzner,  (Ref.  l) .  The 
altitude  h  in  turn  is  defined  by 
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with  re  the  mean  earth  radius ,  and 


r  arf(l*e)/(l  +  ec+4,v). 

The  quantities  Pp  ,  Glp  are  defined  by 

P,  *  r>f  r,1  /i 
&rsrf/>,.  /* 

where  P  is  a  pitching  moment  parameter,  and 
/"*«  is  a  damping  in  pitch  parameter. 

These  parameters  are  defined  and  discussed  in  Volume  I  of  this  series  of  re¬ 
ports.  It  is  desired  to  calculate  9  and  possibly  the  angle  of  attack 
o t  •  0  -  Y  for  "all  values"  of  the  tine  anomaly. 

Even  though  the  differential  equation  is  a  linear  second  order  equation 
it  does  not  lend  itself  readily  to  solution  in  terms  of  elementary  or  tabu¬ 
lated  functions  particularly  since  the  variation  of  density  with  altitude  is 
a  tabulated  function.  Direct  numeric  Integration  appears  to  offer  a  rather 
straightforward  approach,  provided  its  limitation  to  a  finite  range  of  v 
can  be  circumvented. 

The  means  used  to  obtain  results  valid  for  all  vr  constitute  the  first 
section  of  this  report.  The  second  section  describes  the  method  of  numeric 
integration  adopted.  The  final  section  discusses  the  computer  program  written 
to  implement  these  results  while  the  program  Itself  is  given  in  the  appendices. 
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1.  DISCUSSION  OF  THE  LINEAR  SECOND  ORDER  DIFFERENTIAL 
EQUATION  WITH  PERIODIC  COEFFICIENTS 


It  is  the  purpose  of  this  section  of  the  report  to  show  that  three  inte¬ 
grations  carried  out  over  one  orbital  period  suffice  to  furnish  information 
about  the  boundedness  or  stability  of  solutions  of  Eq.  (l) .  Equations  with 
arbitrary  periodic  coefficients  are  discussed  to  insure  independence  of  the 
tabulated  atmospheric  density  function.  In  this  form  the  presentation  is  an 
amplification  of  the  concise  discussion  given  by  Coddington  and  Levinson, 

(Ref.  2)  and  a  minor  extension  as  well  since  the  reference  discusses  homogeneous 
equations  only. 

1.1  General  Case 


Let 


y  +  o Ly  ^  j3> y  -  Y  (5) 

be  a  linear  differential  equation  in  which  c* ,  ft  ,  and  Y  are  periodic  functions 
of  t  ,  with  common  period  of  2rr  .  The  dots  denote  differentiation  with  respect 
to  t  •  Eq.  (l)  is  of  this  type.  Let 

y  +  0 ty  *  y  -  0  (6) 

be  the  associated  homogeneous  equation.  Now  let  y  be  that  solution  of  Eq. 

(5)  that  has  ° 

y^  1 0)  =  0 ,  y  (0)  =  0 

and  let  y  ,  y^  be  solutions  of  Eq.  (6)  that  have  respectively 

yfo)  =/ ,  /  (0)  -  0, 

y  (0)  =  0,  y^  (0)  *  /. 

Since  these  solutions  are  obviously  linearly  independent,  any  solution  of  Eq. 

(5)  may  be  written 

y  •  y„ +  *  b/t  ■  (7) 


For  t*  0,l7T  Eq.  (7)  yields 


I 


f 

I 


y  lo)  -  y  (o)  *  A  y  ( 0)  ■*  t>  y^  10)  r  a 
!/  '  9o  ^  *  &  y  *  i>  y  (o)  8  4  j 


y(2-»r)=  y  (In)  ■*■  0-  y  (In)  +  b  yt  (2  n) 
y  (  tn)  -  yo  (■  2rr  )  ■*■  A  ^  l-2.Tr)  ■+  i>  ^  (  <2ir) 
If  it  is  now  required  that  y  have  period  2  7T  , 

y  (hr)  ~  y  to)  J  y  Utt)  -  y  ( O) 

Eqs.(8)  may  be  combined  to  yield 

y  [In)  =  a.  [l '  y  (2n)l  -  b  y%  Un) 
y  [  2n)  :  “  a. y  Un)  *►  b  [  I  -  yx  (in)] 

Let 

vy  *  f  /,  (l7r>  /,  Un) 

[  y  (In)  yt  (iir) 

then  Eq.  (10)  may  be  written 

y  Un) 


(8) 


( 2rr) 


(9) 


(10) 


A  unique  periodic  solution  to  Eq.  (5)  is  thus  seen  to  exist  is  non- 

singular.  It  will  be  shown  below  that 

cU  (I-  w)  =  0 


implies  a  resonance  in  the  conventional  sense.  In  the  non-resonant  case  let 
the  periodic  particular  solution  be  denoted  by  the  subscript  : 

/f  *  /o  +  + 

where  A  and  b  are  the  solution  to  the  system  of  equations  given  by  Eqs.  (10). 
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Now  define  the  notation 


Then 


y  =  m.AX  (  i/I  )  . 

:  o#t  s  2  rr 

yT  4  y.  *,a,y.  -  IMx , 


(n) 


and  V  exists  and  is  hounded  whenever  y  ,  y  ,  y  ,  exist  on  the  interval 

[0,2irf  ,  and  tit*  H-W)  f  0.  /o  ' 1 

Under  these  conditions,  any  solution  of  Eq.  (5)  may  be  written 


/!/f  +fa.  /,  (12) 

where  a0  and  b0  are  initial  condition  constants.  The  only  change  from  Eq.  (7) 
is  the  choice  of  the  particular  solution.  The  notation 


tn  *  t  -  2h.1T  (13) 

will  be  useful.  In  terms  of  the  conditions  at  t  =  2  h. 7T  Eq.  (12)  becomes 

yM  '  yT  (t n)  +  fl(t 

and  in  turn,  when  t  h  -  2  IT , 

y ( 2w  +1  iff)  *  y  Hit)  +  <*h  y  Hit)  +  bH  yx  ilir) 

s  y?(o)+  y(o)+  b n.,  yx(o)  *y9(o)+  an.,. 

Since  yf  (It r)  =  (0) , 

*  an  y  (Iff)  -  b^yx  l2tr)  • 

Similarly 


V/  *  «»  y.  (2’r)  *bnytUv')- 
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Thus 


•  w  fa.' 

«  w 

IK 

The  solutions  of  Eq.  (5)  remain  bounded  if  W* remains  bounded  for  all  .  This 
in  turn  is  the  case  if  the  magnitudes  of  the  eigenvalues  of W  are  at  most  1. 

In  this  case  an  upper  bound  on  the  amplitude  of  solutions  is 


/  *  7+  +  l7,  * 7 ). 

1.2  Resonant  Case 

If  (I-W)  is  singular,  the  system  of  equations 

O'  (I-w)  \a 

uj 

has  an  infinity  of  non -trivial  solutions.  Let  one  of  these  be 

and  let 

X*  '  a*/.  *  b*Xx  • 

Then  Eq.  (15)  may  he  rewritten 
0  »  (  I  -  w)  [  <k* 
b* 


or 


whence 


b* 


W 


a 

b* 


X  (o)  xl 

‘  A* 

r  lx,  (2*}  Xx{lr'] 

■  A*' 

/,  (°)  Xl  (0,J 

lx,  Uv)  Xl  {l*\ 

b* 

J 

(14) 


(15) 


(16) 


(17) 


7 


so  that 


y  *  (0) 

/H 


y*  Utt)  1 

y*  Uw)J 


(18) 


Thus  y*is  a  non-trivial  solution  of  period  2r  of  the  homogeneous  equation, 
Eq.  (6).  Since  the  period  of  y*  matches  that  of  the  forcing  term  of  Eq.  (5) 
this  meets  the  accepted  notions  of  resonance. 


Conversely,  if  a  solution  of  the  homogeneous  Eq.  (6)  satisfies  Eq.  (18) 
the  argument  may  be  reversed  to  show  that  Eq.  (15)  is  vs.lid. 


In  the  resonant  case,  the  system  of  equations,  Eqs.  (10),  may  or  may  not 
have  a  solution.  If  it  does,  the  analysis  for  the  non -resonant  case  applies, 
except  that  y~  is  no  longer  unique.  If,  however,  the  system  of  equations, 
Eqs.  (10)  is  inconsistent,  no  periodic  solution  of  Eq.  (5)  exists.  In  this 
case  let 


/  =  /o  -  «.  x,  *  b»  /*  U9) 

be  some  solution  of  Eq.  (5),  where  a.  and  k.  are  again  initial  condition  con¬ 
stants.  By  suitable  choice  of  the  initial  condition  constants  Eq.  (19)  may 
be  written 

/  M  =  X.  (tJ  X  «n>  +  /a  LU)  (20) 

where  th  is  defined  in  Eq.  (13 )  •  Now 

y  Un  *  2  Kir)  =  yo  f In  )  +  y  (in  )  +  ^  Un) 

=  /•  {o)  +  X  l0)  *  Xa  l0)  s  «*♦,  , 

and 

y  (2n  -  2h  rr)  *  ^  (2v)  ♦  lilt)  »  £  Uwl 

‘  y.  i°>  *  «...  /,  10)  -  io> .  . 

Thus 


8 


and,  proceeding  recursively 

s  f/w  ^7r;l  +  w  f  [ /t  Uir)  1  +  W  '  au\ 

.  b-J  [K  Urr)l  [l/.(2ir)J  *>-]) 

s  (1  +W)  f  (2ir)  +  Wl 

l  K  <**> 

»•••-*  (  £  W  * )  f  y  (nr)  +  W  **' 

i*o 

K  lln{ 

Now  the  vector  [ak+i  L*,,  J  can  remain  bounded  as  *v  increases  only  if  the  sum 

5n  *  H  W*  (22) 

t*o 

remains  bounded. 

Let  <rnj  be  the  eigenvalues  of  .S*  and  u /j  be  those  of  W ,  where  t  *  /,  1  . 
The  eigenvalues  of  the  metric  polynomial,  Eq.  (22)  have  the  relation 

rKl  t  fi-  w;')/(i-w;)  ,  w;  *  o,i 

~  0  ,  Wi  -  0 

»  n  )  **1  *  i 

Thus,  the  definition 

W  *  hVfbK  (  I  W,  I  j  I  w,  I  )  (23) 

and  relation 


0  <  v  <  / 


(24) 
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imply 

I  P'n'i  I  *  I  J/  (  /”  )  I  *  I  /  (  I  ~  Ur )  aa.  —+  oo  . 


But  if  the  relation,  Eq.  (24),  is  satisfied  the  second  term  of  Eq.  (2l) 
is  also  bounded,  and  the  amplitude  of y  is  limited  to 

I y\  i  L  ^ y,  tin)  1 ~+  y  l  Lir)^  *  ly  +  y )  V  *.1  +  bj^  (25) 

I  -  vr  *  * 

where  the  bars  are  defined  in  Eq.  (ll). 

The  derivation  of  the  bound,  given  by  Eq.  (25),  did  not  use  the  fact  of 
resonance,  thus  it  is  applicable  even  to  the  non-resonant  case,  although  in 
the  latter  case  the  bound  given  by  Eq.  (l4)  will  usually  be  sharper,  (it 
should  be  borne  in  mind  that  the  symbols  a0  and  have  different  significances 
in  the  two  expressions.) 


1.3  Stability  Criteria 

The  eigenvalues  of  the  matrix  W  satisfy  the  equation 


0  ■ 


’  -  yx  Uv)  -/x 
<yt  Utt)  w  *  y  Uv) 


u?1  -  f  yt  (in')  ♦  y  ( 2ir) J  w  ♦  W 


so  that 


w  s  y  fin)  ♦  y  ( Irt) 
~~  ~ 


/,  (_^lz  A  Un) 


(26) 


When  the  eigenvalues  are  complex,  they  are  conjugates,  so 

i* 


3 


y  (2ir)  4  y  (2rr) 
l 

^  W  , 


+  M  w  - 


j  y  Uv)  -*  y  Uw) 

L  ~  1  . 


and  thus 

I  w,  I  * 


1 


When  the  eigenvalues  are  real,  let  the  one  larger  in  absolute  value  be  . 

Then 

I  u>,  I  l  y  ( 2 tr )  *  yt  (In)  i  */  eh*  W*  *  - 

l 

For  bounded  motions,  therefore,  it  is  necessary  that 

cUi  w  *  l. 

It  is  to  be  noted  that  tUt  W  is  the  Wronskian  of  the  system  ( y  }  y  )  , 
evaluated  at  t  *  2  TT  .  Abel '  s  formula  for  the  Wronskian  is :  1 

UTu)  *  WLO)  e  ult)clt , 

Thus 

^eJtW  =  W"( 2  n)  c  yt  (o)  yt(o) 

y.  <°>  h  <°> 

Now  if  oi  is  identically  zero,  or  even  if  it  has  an  average  value  of  zero  on 
the  interval  (0,Ztr), 

clitw  *  /. 

Then,  referring  to  Eq.  (23), 
vr  >  /. 

It  may  be  concluded,  therefore,  that  in  the  absence  of  the  damping  term,  the 
resonant  solutions  are  unbounded  unless  Eqs.  (lO)  happen  to  be  dependent.  The 
non-resonant  solutions  were  shown  above,  Eq.  (14),  to  be  bounded,  unless  yjr  >  /. 

1.4  Short  Term  Boundedness 

In  the  application  of  the  results  obtained  above  to  Eq.  (l)  it  must  be 
recalled  from  Volume  I,  Appendix  C  of  this  series  of  reports,  that  the  con¬ 
stants  (or  elements)  describing  the  orbit  are  themselves  subject  to  drift  as 
the  orbit  decays .  Thus  the  results  are  applicable  to  some  large  but  finite 
number  of  periods.  Under  such  a  limitation  an  ultimately  unbounded  motion 


e  /.2Tr*(t)</t  ,  e/0  ollt)  cit  m 
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may  still  be  acceptable  if  its  rate  of  growth  is  small.  For  the  nth  period 
the  bound  expressed  by  Eq.  (25)  may  be  written: 

y  *  urK-  1  V yc  ( lit)1  *  ■*  w  **'  ly  *  y  )  y  (o)v+j/(o)  |27 ) 

us  -  I 

while  in  the  non -resonant  case  the  bound  given  by  Eq.  (l4)  becomes 

y  (tK_' )  sy  +  w  K*'  J  i  y^io)-  y(o)]1  *  l  yr  (o)  ~ylo)]1  .  (28) 

Now  the  total  lifetime  of  the  satellite  may  be  divided  into  regions  during 
each  of  which  the  parameters  of  the  orbit  may  be  considered  constant,  and 
during  each  of  which  the  growth  or  decay  of  the  pitch  amplitude  may  be 
estimated  by  means  of  Eq.  (27),  or  in  a  non -resonant  case  by  Eq.  (28). 

I.5. A  "Built-In"  Check  on  the  Accuracy  of  the  Numeric  Integration 

Let  Eq.  (5)  have  a  unique  periodic  solution  y’  and  let 

/  T 

yr  *  ye  *  /, 

where  yg  is  the  even  part  of  y_  and  y*  is  the  odd  part  of  y ^  .  Further,  sup¬ 
pose  that  in  Eq.  (5)  the  functions  ol  and  V  are  odd  functions  of  t  ,  and  /ft 
is  an  even  function  of  t  .  Then  Eq.  ( 5 )  becomes 

( /e  *  *  /« *  x.  * + '  >v  * *  y- 

The  parity  properties  of  functions  guarantee  that  the  first  parenthesis 
on  the  left  is  an  even  function  of  t  ,  while  the  second  parenthesis  is  an 
odd  function.  Since  the  right  member  is  an  odd  function,  the  first  parenthesis 
must  vanish.  If  there  existed  a  non-trivial  even  function  y  the  uniqueness 
of  yr  would  be  destroyed.  Thus  * 

y,  *°> 

and  y_  is  an  odd  function,  so  that  y  lo)  ~  0  .  Then  in  the  solution  of  Eqs. 

(10I : 


a  -*  0 


yt  ■  /.  -  • 

This  situation  arises  in  Eq.  (l)  when  the  damping  in  pitch  parameter, 
r%  ,  vanishes.  In  this  case  the  size  of  4  computed  from  Eqs.  (10)  becomes  a 
measure  of  the  accuracy  of  the  numeric  integration  used  to  obtain  the  coeffi¬ 
cients  in  these  equations. 
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2.  INTEGRATION  PROCEDURE 


The  choice  of  integration  procedure  involves  a  balance  between  accuracy, 
speed,  stability,  and  auxiliary  "housekeeping"  requirements.  Gill,  (Ref.  3) 
has  derived  a  very  elegant  version  of  the  Runge-Kutta  fourth  order  procedure 
in  which  the  housekeeping  is  minimized,  in  that  only  one  storage  cell  need 
be  assigned  to  each  variable  to  be  integrated.  Milne  (Ref.  4,  Section  38) 
earlier  gave  a  predictor-corrector  method  which  obtains  the  same  accuracy  with 
half  as  many  substitutions  into  the  differential  equation,  but  requires  that 
a  history  of  past  values  of  the  variables  be  kept.  Milne's  predictor  formula 
requires  a  history  of  four  previous  values  of  the  derivative,  and  if  it  is 
desired  to  double  the  step-size  in  the  course  of  the  integration  a  history  of 
seven  previous  values  must  be  available,  along  with  a  record  of  the  time  at 
which  the  step  size  was  last  doubled.  The  starting  procedure  for  the  predictor- 
corrector  method  is  necessarily  relatively  involved.  Since  the  corrector 
formula,  Simpson's  Rule,  requires  a  history  of  only  two  previous  values,  it 
would  be  desirable  to  have  a  predictor  formula  with  an  equally  short  time-base, 
even  at  the  cost  of  some  accuracy  and  stability.  Milne  (Ref.  4,  Sections  30 
and  3l)  gives  a  procedure  for  deriving  such  formulas  to  arbitrary  specifications 
and  for  estimating  their  accuracy.  The  procedure  is  here  used  to  obtain  the 
desired  formula. 

2.1  Derivation  of  Predictor  Formula 


In  the  identity: 


y  (t0 +1)  -  y  (t,-£)  =  ^  y  dt 


(29) 


the  assumption  that  the  second  term  on  the  right  is  a  linear  combination  of 
y(t0),  y(t'-Jl),  A  yttm)  and  4  y  ( t,  -  J>  )  ,  may  be  written 

I  y  Jit  *  A  y  (t9)  +  b  y(t0~A)  +  C  A  y  (t9)  +  D  A  y  (t9~  A)  .  (30) 

"  -4 

The  coefficients  A,  B,  C,  D  may  be  determined  by  requiring  that  Eq.  (30)  be 
exact  for  the  first  four  terms  of  the  Taylor  expansion  of  about  .  Ihus, 
let 


y  (t) 

where 


«■  a ,  (t  -t 0)  «■  at  It' t.)1 

Z 


a5  (t-t.) J  «•  - 
6 


(31) 


<JK y  (t.) 


^  * 
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Substitution  of  Eq.  (3l)  into  Eq.  (29)  and  Eq.  (30)  gives  respectively 

y  Jt  s  l  a.,  £  +  j 


1 


t 

to  +1 


J*  y  it  «  A«.  *8  la.  -  aj  *  a,  ll  -  *>£') 


z 

Equating  the  coefficients  of  ao>& iAx>  and  a»  respectively,  yields: 

A+  8  =  0  1 

C  -  D-  8  =  2 

ID-J6=  J 

The  solution  to  this  system  is 

A  =  ~  ‘t  j  6  *  C- V,  D  3  2  . 

Then  Eq.  (30)  becomes 


I 


y  Jt  a  -h-ylt.)*'*ylt.-4)  *  'fiy(t.)*lf,yit.-Jt), 

t ,  *  J 

which  yields  in  Eq.  (29): 

y(t.^)  *  y(t,-i)*h-{y(t„~i)-y(t.l  +  d[y  It.-i  I  *2  y(t,)]} 

3  /,  (t.  *4 ) 


(32) 


where  the  subscript  p  stands  for  "predicted  value".  Its  derivation  insures 
that  Eq.  (32)  is  exact  to  third  order  terms. 


Ohe  Taylor  series,  Eq.  (3l)  may  be  made  exact  by  including  the  residue 


term: 


y  (t)  -  a,  +  a,  (t- 1.)  +  a ,  (t-tj*  *  at  (t-t.)}  +  .  iLllsi!  (33) 

l  6  clt+  2V- 
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where  s  is  on  the  interval  L  t0j  t  ]  . 

Now  the  error  committed  "by  the  use  of  yt  (t0  +  ^>)  in  place  of  y 
may  he  obtained  by  subtracting  the  right  member  of  Eq.  (32)  from  its  left 
member.  If  each  term  in  this  difference  is  expressed  in  terns  of  Eq.  (33) 
or  its  first  derivative,  the  coefficients  of  a,,  ,  A*  ,  and  A3  will  be 

zero  because  Eq.  (32)  is  accurate  to  third  order  terms.  The  residue  terms 
will  remain: 

+  d\u.)  -5A*ytsx)  +  8  A*yts1t]  , 

1*1  At*  At*  It*  J 

where  ,s(  is  between  t0  and  t8+^,  while  and  are  between  tp  and 
The  three  distinct  values  arise  from  the  fact  that  f  in  Taylor's  Theorem, 

E<1-  (33),  cannot  be  expected  to  be  the  same  for  the  expansions  of  yit0  +  A), 
yU9-  A),  and  y  lt0- /) .  However,  Eq.  (32)  may  be  shown  to  satisfy  Milne's 
Theorem  2,  (Ref.  k,  Section  3l)  which  in  effect  permits  the  replacement  of  S,  , 
S j  ,  and  by  a  single  5  between  tt-Jl  and  The  error  becomes 

A*yis>  V  , 

7F - T~ 

and  Eq.  (32)  may  be  rewritten: 

y(t*  *S\)  *  y(t0-^)  +  4-{y(t0-i) -y  (te) 

+  A[y(t0-£)  +  2y(t.)]}  +  d+y  . 

I  IF'  6 

2.2  Corrector  Formula. 

Simpson's  rule  is 

ylt.*t)*  ylt.-J!)+  Uytt.-&)+*y(t.)+ylt.+M  --  yftt.*A) 

3 

where  the  subscript  3  denotes  "calculated  from  Simpson's  rule".  An  error 
analysis  for  Simpson's  rule  yields 


(3fc) 


(35) 


y  (t,  *1)  *  y  (te  -  A)  +  &  [y  *  ky  (tp)  +y(t0+Ji)]-  JSy  .  (36) 

•  3  dts  ’  90 

In  the  derivation  of  the  error  term  in  Eq.  (35)  the  assumption  was  made  that 
the  other  terms  on  the  right-hand  side  are  known  precisely.  If  this  assumption 
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is  replaced  by  the  assumption  that  y{ to+^Jwas  itself  computed  by  Eq.  (32; 
the  error  term  is  larger,  in  fact 

y(to+I)=y(to-4)  +  l[y(t0-t)^y(t.)+yTlt.+Ml  +  l*-&  If  .  (37) 

3  Jit5  *  90 

If,  however,  y^  (t0+A)  is  used  in  place  of  y  (t0+i>)  the  error  term  is  about 

6  d5y  . 

5  *  3t^“? 0 

2.3  Integration  Scheme 


Kie  predictor- corrector  scheme  is  then: 


a)  With  y  *  0*  use  Eq.  (3*0  to  predict  0'(v  +  Av)  from  0'lv),  0’(v-Air), 
@’*( v),  &”lv-Av). 


b)  Withy*  0  use  Eq.  (36)  to  predict  0(v*Av)  from  0  (if- A  if), 
0’(v-av)  ,  0’ (v)  and  the  predicted  9'(v+av). 

c)  Use  the  predicted  fl(v  +  av)  and  0’tv^Av)  in  Eq.  (l)  to  predict 

e"lv  +  AV)  . 

d)  Use  the  predicted  0"(v+Av)  in  Eq.  (35)  with  y*0’  to  calculate 
the  corrected  0'(v  +  Av). 

e)  Use  the  corrected  ©‘(v  +  Av)  in  Eq.  (35)  withy  e  0  to  calculate 
the  corrected  0  ( v  A  v) . 


f)  Use  the  corrected  ©(v+axr) 
the  corrected  0"(v-*-av). 


and  &’  (\r  +  A  vr) 


in  Eq.  (l)  to  calculate 


The  small  size  of  the  error  term  in  Eq.  (36)  makes  it  unlikely  that  an 
iteration  of  steps  d,  e,  f  would  improve  the  result,  if  the  prediction  of  steps 
a,  b,  c  is  at  all  adequate.  This  adequacy  may  be  tested  by  comparing  the 
predicted  and  corrected  values  of  0  ”  ( v  a  v)  .  Since  the  error  estimates  for 
the  prediction  exceed  those  for  the  correction  the  difference  between  predicted 
and  corrected  values  is  a  pessimistic  estimate  of  the  error  incurred  at  each 
step  of  the  integration.  Let  t  (av)  be  the  magnitude  of  the  difference  between 
the  predicted  and  corrected  values  of  9*  (v  +  tv)  ,  and  let  6  be  the  maximum 
tolerable  difference.  Then  if 


£U\r)  *  £/3Z 


the  error  estimates  for  Eq.  (35)  suggest  that 


£  (1  AV)  *  €  • 


1  / 
i.O 


Thus  a  stiff iciently  small  £  justifies  a  doubling  of  the  step  size.  If,  how¬ 
ever 


£  (  a  \r)  >  £. 

the  integration  must  be  restarted  with  a  smaller  value  of  iV  . 

Starting,  as  well  as  restarting  may  be  accomplished  by  Taylor's  formulas 

6'(v+&yr )~  0’(v)  +  ©"(vjAV"  +  d39(s)  (av)1  ,  (38) 

dv>  *  2 


9(v  +  av)  -  9lv)  *  @'(v)av  +■  9n(yr)(A\r)^  +  di  Bis)  ( Air)^  ,  (39) 

2  dv*  6 

where  j  is  between  \r  and  ir  +  AV  .  The  relatively  low  powers  of  AV  appear¬ 
ing  in  the  error  terms  of  Eq.  (38)  and  (39)  suggest  that  small  values  of  aw 
be  used  to  retain  accuracy.  Since  the  predictor-corrector  scheme  makes  it 
easy  to  double  the  step  size  there  is  but  small  penalty  for  setting  the  initial 
value  of  aw  too  small. 

Computational  experiments  have  been  conducted  comparing  the  procedure 
presented  above  with  Gill's  version  (Ref.  3)  of  the  Runge-Kutta  fourth  order 
procedure  and  with  the  Kutta  fourth  order  formula  as  quoted  by  Milne  (Ref.  5)* 
All  three  processes  gave  essentially  the  same  accuracy  for  the  same  step  size; 
as  mentioned  above,  the  present  process  requires  half  as  many  substitutions 
into  the  differential  equation.  It  should,  however,  be  noted  that  the  rela¬ 
tively  large  coefficients  appearing  in  Eq.  (3^)  tend  to  favor  the  build  up 
of  "inherited”  errors,  so  that  the  present  method  cannot  be  expected  to  be  as 
stable  as  the  Runge-Kutta  process  for  differential  equations  with  unstable 
solutions . 

The  potential  build-up  of  inherited  error  is  not  of  concern  in  the  present 
study  for  three  reasons .  Firstly,  the  range  of  integration  is  relatively 
short  so  that  inherited  error  cannot  build  up  very  far.  Secondly,  the  most 
critical  cases  are  the  most  unstable  cases,  but  in  these  the  demonstration 
of  instability  is  obvious.  Thirdly,  application  of  the  accuracy  criterion 
developed  in  Section  1.5  to  the  computer  results  justifies,  a  posteriori,  the 
assertion  that  the  error  build-up  was  nil  to  four  decimal  places. 
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3-  COMPUTER  PROGRAM 


Program  Descriptions  of  the  various  components  of  the  Computer  Program 
written  to  calculate  solutions  of  Eq.  (l)  are  included  in  Appendices  A  to  F. 
Sample  input  and  output  data  are  in  Appendix  G.  The  program  has  been  checked 
out  and  used  at  the  North  American  Aviation,  Inc.  Columbus  IEM  JOk  computing 
installation  to  generate  the  numerical  results  discussed  in  Volume  I  of  this 
series  of  reports.  Certain  systems  features  of  this  program,  such  as  the 
choice  of  input  medium  (magnetic  tape  or  punch-cards),  assignment  of  magnetic 
tape  unit  numbers,  labelling  of  library  routines,  etc.  may  have  to  be  varied 
if  this  program  is  to  operate  successfully  at  other  installations  or  on  other 
computing  machines.  It  is  believed  that  the  information  furnished  is  adequate 
to  enable  anyone  familiar  with  the  IEM  Fortran  Manuals  (for  instance  Ref.  6 
and  7) ,  and  the  operating  system  of  his  installation  to  adapt  the  source 
language  deck  to  that  installation. 

The  program  was  written  in  two  portions.  The  first  portion  integrates 
the  differential  equation,  Eq.  (l),  and  its  homogeneous  part,  and  records  the 
results  on  a  magnetic  tape  as  well  as  on  the  output  medium.  The  second  portion 
reads  the  magnetic  tape  and  produces  punched  cards  suitable  for  automatic 
plotting.  The  portions  are  discussed  separately  below  and  a  brief  discussion 
of  the  atmosphere  table  look-up  routine  is  also  given.  This  section  concludes 
with  some  remarks  on  a  potential  application  of  the  computer  program. 

3-1  Integration  Portion  of  the  Program 

All  input  data  is  read  by  means  of  a  "sparse  data"  routine  so  that 
information  unchanged  from  one  case  to  the  next  need  not  be  written  or  key¬ 
punched  repetitively.  The  only  item  of  input  data  modified  by  the  program 
is  the  identification  number  used  to  identify  records  on  the  intermediate 
storage  tape,  which  is  automatically  advanced  from  record  to  record.  To 
further  facilitate  the  preparation  of  input  data  the  primary  input  variables 
are  presented  in  lists;  the  program  systematically  uses  each  combination  of 
list  entries.  All  angles  in  input  and  output  are  expressed  in  degrees; 
however,  inside  the  computer  and  on  the  intermediate  storage  tape  they  are 
carried  in  radian  units. 

In  the  integration  of  the  differential  equations  care  is  taken  to  provide 
common  abscissas  from  one  integration  to  the  next.  This  is  accomplished  by 
using  a  computing  interval  obtained  from  the  print  interval  by  successive 
halving,  and  by  not  doubling  the  step  size  (when  doubling  is  permitted  by  the 
criteria  of  Section  2)  if  this  would  cause  the  end  of  the  print  interval  to 
fall  inside  a  computing  interval.  The  accuracy  of  the  integration  is  further 
insured  by  a  test  which  aborts  the  calculation  if  the  step  size  does  not 
expand  in  the  first  three  intervals.  If  the  step  size  fails  to  expand,  the 
accuracy  of  the  starting  formulas  Eqs.  ( 38) ,  (39)  is  less  than  that  implied 
by  the  predictor- corrector  method.  The  calculation  should  be  attempted  again 
with  a  smaller  minimum  step  size.  The  integration  is  also  aborted  if  the 
magnitude  of  any  of  the  0's  exceeds  10?,  such  an  event  would  imply  a  gross 
data  error,  or  a  machine  failure.  The  number  of  computing  steps  taken  during 


each  print  interval  is  printed  out  so  that  the  step  size  used  may  be  deduced, 
the  total  number  of  steps  for  each  integration  is  printed  out  as  an  aid  to 
future  computing  time  estimates.  As  compiled  for  the  IBM  70b,  the  program 
will  integrate  about  8.5  steps  per  second. 

3.2  Card -Punching  Portion  of  the  Program 

This  section  may  be  combined  with  the  first  section  in  a  single  computer 
run,  or  it  may  be  used  separately  if  the  intermediate  storage  tape  is  saved 
from  run  to  run.  The  second  option  permits  inspection  and  editing  of  the 
curves  to  be  punched  into  cards.  Since  the  on-line  punching  of  cards  is  a 
relatively  expensive  process  the  card  format  was  designed  to  admit  up  to 
fifteen  curves  per  card  deck,  while  retaining  complete  flexibility  as  to  size 
and  placement  of  the  plots.  This  portion  also  reads  its  input  data  (specify¬ 
ing  the  editing  information)  with  the  "sparse  data"  routine.  The  only  input 
item  modified  by  the  program  is  the  card  deck  identification  number  which 
serves  as  a  signal  that  all  information  to  be  punched  has  been  accumulated. 

This  number  is  reset  to  zero  after  the  card  deck  has  been  punched  to  indicate 
that  the  store  of  information  is  again  incomplete. 

The  punched-card  output  for  automatic  plotting  assumes  a  plotting  device 
capable  of  accepting  one  coordinate  pair  per  card.  Abscissa  and  ordinate 
may  be  placed  in  any  of  sixteen  four  digit  fields,  or  the  abscissa  may  be 
suppressed  in  favor  of  an  abscissa  punched  for  another  curve.  The  abscissa 
field  will  contain  the  independent  variable  v'  .  The  ordinate  may  contain 
any  function  of  the  form 

J-  *  I  ur^l  0,  -*■  ur,  0,  t  Bj,  *  ury  Y 

where  &0  is  the  particular  solution  of  Eq.  (l)  corresponding  to  y  of 
Section  1, 

0i,  are  the  solutions  of  the  homogeneous  part  of  Eq.  (2)  correspond¬ 
ing  to  y{  >  yx  °*  Section  1, 

Y  is  defined  in  Eq.  (2) 
and  the  u^'s  are  input  data. 

For  vr  *  0 ,  urT  »  /  ,  f  will  be  a  solution  of  the  differential  equation,  Eq.  (l). 
When  wy  •  -1 ,  /  will  be  the  angle  of  attack,  «<  .  When  u =  - 1,  the  input 
values  of  u t,  and  u^will  be  ignored,  instead  the  program  will  use  the  pre¬ 
computed  values  needed  to  make  f  periodic,  of  orbital  period. 


The  card  decks  produced  will  contain,  in  addition  to  the  abscissa -ordinate 
pairs,  four  terminating  cards.  The  first  two  of  these  are  calibration  check 
cards  which  will  position  the  plotting  head  to  the  lower  left  hand  and  upper 
right  hand  corners  of  the  plot  respectively,  the  third  contains  an  11 -punch  in 
Column  1  which  may  be  used  to  stop  the  plotting  device,  and  the  fourth  is 
completely  blank,  serving  as  a  separator  between  decks.  The  body  of  the  deck, 
and  the  first  two  terminating  cards  have  a  O-punch  in  Column  1,  sequence  number 
in  Columns  5*8>  and  a  deck  identification  number  in  Columns  2-b.  The  four 
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digit  fields  consisting  of  Columns  9-12,  13-16,  etc  to  69-72  are  numbered 
F  =  1,  2,  l6  respectively  and  are  the  data  fields  proper  (Columns  73-80 

are  not  readily  accessible  to  the  IBM  JOk  computer) . 

The  tape  file  of  data  is  searched  completely  for  the  record  required  by 
the  input  data,  so  that  in  principle  these  records  could  be  demanded  in  any 
order.  The  search  is  most  efficient  however  if  the  records  are  demanded  in 
the  order  in  which  they  are  stored  on  the  tape. 


Atmospheric  Properties  Subroutine 


The  atmospheric  properties  routine  used  in  the  program  is  capable  of 
simulating  atmosphere  tables  based  on  the  ideal  gas  law  and  on  linear  varia¬ 
tions  of  temperature,  gas  constant,  and  specific  heat  ratio  with  geopotential 
altitude.  In  the  present  program  the  routine  was  used  to  obtain  tne  density 
distribution  of  the  1959  AjRDC  (Minzner)  Atmosphere,  (Ref.  l),  which  was 
represented  by  use  of  the  Molecular  Scale  Temperature  in  place  of  the  usual 
temperature.  The  detail  input  data  are  given  in  Appendix  G. 


In  its  present  form  the  program  can  furnish  information  about  the 
stability  in  pitch  of  a  given  satellite  in  a  given  orbit.  If  the  detailed 
motion  is  not  of  interest,  only  the  first  portion  of  the  program  is  required. 
For  a  design  study,  .knowledge  of  the  stability  boundaries  in  various  parameter 
planes  of  the  differential  equation,  Eq.  (l),  would  be  desirable.  These 
occur  when  the  magnitude  of  the  larger  eigenvalue  ui  given  by  Eq.  (26)  passes 
through  the  value  1.  With  minor  modifications  the  subroutine  INTEG  of  the 
first  part  of  the  present  computer  program  could  be  used  as  the  basis  of  a 
program  to  search  iteratively  for  these  stability  boundaries.  Points  of 
departure  for  these  searches  could  be  obtained  from  the  eigenvalues  of  the 
sine  problem  to  which  the  present  problem  reduces  itself  when  the  orbital 
eccentricity  is  zero. 


'—.jawtuvnai  ....  -wta-tSRsi 
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A  NOTE  ON  THE  APPENDICES 


Appendices  A  to  E  are  presented  in  the  format  of  a  standard  SHARE  program 
description . 
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APPENDIX  A 


PROGRAM  DESCRIPTION  -  PITCHING  MOTION  OF  A  SATELLITE 
NORIK  AMERICAN  AVIATION,  INC.  COLUMBUS 
ENGINEERING  PROGRAM  DESCRIPTION 


1.  Identification 

a.  Pitching  Motion  of  a  Satellite,  IM  091 

b.  0.  C.  Juelich,  December  i960 

c.  Research  Group  -  350 

d.  Fortran  Source  Deck  is  up  to  date. 

2.  Purpose 

Numerically  integrates  one  particular  and  two  homogeneous  solutions 
of  the  pitch- equation.  Solutions  are  printed,  and  optionally  punched  into 
cards  in  scaled  form  suitable  for  automatic  plotting. 

3-  Restrictions 

a.  Uses  Tape  5  as  input  tape 

Tape  6  as  printer  output  tape 

Tape  2  for  intermediate  storage  of  results 

Tape  7  to  store  part  of  the  program 

b.  Operates  in  N.A.A.  monitor  for  IBM  704  Fortran 

4.  Method 

a.  The  mathematical  method  is  described  in  Report  WADD  TR  61-133  Vol. 

Ill,  of  which  this  description  is  a  part. 

b.  The  program  is  designed  to  terminate  calculations  if  the  step  size 
has  not  adequately  expanded  in  the  first  three  print  intervals  or  if 
any  of  the  solutions  exceed  10  million  degrees  in  amplitude . 

c.  Punched  Card  output  is  described  in  Section  3*2  of  this  report. 

5.  Use 


The  program  consists  of  two  parts.  The  first  part  generates  solutions 
of  the  differential  equation  and  records  these  on  magnetic  tape.  The  second 
part  retrieves  selected  solutions  and  scales  them  for  automatic  plotting. 

The  two  parts  can  be  combined  in  one  machine  run  or  they  can  be  used  on 
separate  occasions  if  the  intermediate  storage  tape  (tape  2)  is  preserved. 
The  latter  alternative  permits  inspection  of  the  printed  results  before 
card- punching  is  undertaken. 

Coding  Information 


6. 


a.  The  first  part  consists  of  a  main  program  and  two  subroutines  labeled 
INTEG,  DIFFC.  In  addition  subroutine  AT*0S  (1F113)  and  Utility  File 
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subroutine  FXRCH  are  used,  as  well  as  Library  subroutines  SQRT,  CfiS, 
SIN,  CHAIN,  and  Fortran  system  routines. 

b.  The  second  part  is  Chain  1  on  Tape  J.  It  consists  of  a  main  program 
and  a  subroutine  labelled  SCALE.  In  addition  Utility  File  Subroutines 
FXRCH,  NPSCL,  P/RITE,  are  used,  along  with  Fortran  system  routines. 

c.  Descriptions  of  subroutines  AT!$S,  FXRCH,  NPSCL,  and  WRITE  are  given 
separately. 


Input  Data 


a.  All  input  data  to  both  parts  of  the  program  is  read  by  the  Uility 
File  Routine  FXRCH  into  the  Fortran  C0MM$N  Data  Region. 


b.  Part  one  requires  the  following  data: 


Location  Symbol 
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ration,  Description 


Minimum  computing  step  size  bound,  degrees 
The  minimum  computing  step  size  is  the  largest 
number  of  the  form  d  /lK  (  »\  integral) 

and  &  dv*.  r 


Print  Interval,  degrees 
Should  be  a  factor  of  36o°- 


28  RECID  Record  Identification  Number  or  Control 

Number . 

If  this  number  is  negative  or  zero  the  remain¬ 
ing  data  set  is  ignored.  If  RECID  is  nega¬ 
tive  part  II  is  loaded  and  given  control. 

If  RECID  is  zero  tape  2  is  rewound  and  com¬ 
putation  ends. 

If  RECID  is  positive  it  becomes  the  identifica¬ 
tion  number  of  the  data  case,  if  it  is  omitted 
the  value  of  RECID  from  the  previous  data  case 
increased  by  1.  is  used. 


30  tf  Radius  of  Earth,  feet 

Used  in  Eq.  (4),  but  not  in  subroutine  AIM0S. 

31  K  Product  of  Universal  gravitational  constant 

and  mass  of  earth,  ft3/sec2 
Used  in  Eq.  (3). 

32  TEL  Table  used  for  routine  A2*fte  (1F113) 

See  Appendix  B. 
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Location 

Symbol 

Designation,  Description 

201 

Number  of  orbit  inclinations, 

202 

i  i 

first  orbit  inclination,  degrees 

203 

‘a 

second  orbit  inclination 

et  seq. 

211 

Number  of  perigee  altitudes, 

212 

first  perigee  altitude,  feet 

213 

second  perigee  altitude 

et  seq. 

221 

Number  of  eccentricities, 

222 

first  eccentricity,  dimensionless 

223 

second  eccentricity 

et  seq. 

231 

Number  of  pitching  moment  parameters, 

232 

rt 

first  pitching  moment  parameter 

233 

f'l 

second  pitching  moment  parameter 

et  seq. 

24l 

Number  of  damping  in  pitch  parameters 

242 

rV 

first  damping  in  pitch  parameter 

243 

ru 

second  damping  in  pitch  parameter 

et  seq. 

251 

Number  of  inertia  ratio  parameters 

252 

M, 

first  inertia  ratio  parameter 

253 

M, 

second  inertia  ratio  parameter 

et  seq. 

The  hi.  h*  n«  Mr  nr>-,vM  combinations  of  the 
listed  values  of  i,  h,  e,  0  »  M 

will  be  used  systematically.  The  results 
will  be  identified  by  successive  values  of 
RECID. 


c.  Part  two  requires  the  following  data: 

Location  Symbol  Designation,  Description 


1. 

RECID 

Identification  number  of  tape  record  to  be 
plotted. 

2. 

Weight  of  function  &0  in  curve  to  be  plotted 
Should  be  1,  0,  or  -1. 

3- 

U ry 

Weight  of  V ixr)  in  curve  to  be  plotted 

Should  be  0  or  -1. 

4. 

Weight  of  function  9 ,  in  curve  to  be  plotted 
Ignored  when  wy  =  -1. 

5- 

Weight  of  function  9t  in  curve  to  be  plotted 

Ignored  when  w y  =  -1. 
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Location  Symbol 

6.  /  Imm< 


7*  C. 


8. 


f 


>n«M 


9- 


10. 


c 


11. 


F/ 


12. 


15. 

16. 


* 


17. 


18. 


OIL 


Designation,  Description 

Minimum  value  of  function  in  plotting 
interval,  degrees 

If  f  is  less  than  this  value  the  point  will 
be  plotted  on  the  bottom  margin  of  the  paper. 

Plotter  counts  for  f  . 

Maximum  value  of  function  in  plotting 
interval,  degrees 

If  f  exceeds  this  value  the  point  will  be 
plotted  on  the  top  margin  of  the  paper. 

Plotter  counts  for  . 

Plotter  counts  for  top  margin  of  the  paper. 
The  bottom  margin  is  at  zero  counts. 

Number  of  the  card  field  in  which  f  is  to 
appear. 

Should  be  between  1.  and  l6. 

Minimum  value  of  vr  in  plotting  interval, 
degrees 

If  \r  is  less  than  this  value  the  point  will 
be  plotted  on  the  left  margin  of  the  paper. 

Plotter  counts  for 

Maximum  value  of  v  in  plotting  interval, 
degrees . 

If  v  exceeds  this  value  the  point  will  be 
plotted  on  the  right  margin  of  the  paper. 

Plotter  counts  for 

Plotter  counts  for  right  margin  of  the 
paper. 

The  left  margin  is  at  zero  counts. 

Number  of  the  card  field  in  which  \r  is  to 
appear,  or  control  number. 

If  Fv  is  zero,  V  will  not  be  entered  on 
the  punched  cards. 

Identification  number  for  card  deck  to  be 
punched,  or  control  number. 

If  CID  is  zero  the  scaled  data  are  accumu¬ 
lated  in  core  storage  for  later  punching. 


<WH 


Location 


Symbol 


If  CID  is  not  zero  the  accumulated  information 
is  punched  into  cards  and  CID  becomes  an  ident¬ 
ification  number  in  columns  2  to  k  of  each 
card.  After  punching  CID  is  reset  to  zero. 


d.  Sample  Input  Data  with  resulting  output  follow  the  compiled  listings  of 
the  components  of  this  program. 


Note:  The  program  symbol  for  v  is  PHI,  for  thX  (I"W)is  DEI/TA. 


Main  Program  for  Part  I 
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s  aw-B  *>  .anwa&euns 


•e<Ws»  hk<  *eii«i.  uasMiaima 


Subroutine  INTEG 


Write  Parameter  Combination 
Step  Size 
RECID 

on  Ojtput  Tape 


Coll  DIFFC  to  obtain 
Coefficients  of 
Differential  Equation, 
Calculate  2nd  Derivatives 


Set  working  step  size  to  minimum 
Set  indicator  to  "Retain  step  size" 
Use  Taylor's  Formulas  to  (re)start 
integration 


Convert  Print  Step 

Minimum  Computing  Step 
to  radians 


Call  DIFFC  to  obtain  coefficients 
of  Differential  Equation  at  "next 
point".  Calculate  2nd  derivatives 


Redetermine  Minimum  Computing 
Step  as  (Print  Step)  /2n 


Set  up  Constants  to  test  for: 
Time  to  Print,  Time  to  end 
integration 


Set  up  Test  numbers  for 
Blocks  53A  and  54 


Initialize  Step  Counters 
Print  Counter 


Calculate  Perigee  values  of 
r,  v,  o,  Q,  P,  Vfr/V; 

Set  initial  values  for  integration 


Examine  Step  Size  Indicator 


Save  largest  values  of  0'o 
as  9  max,  0‘  max 


Are  0  max,  0'  max  excessive? 


Advance  Step  Counters 


Have  three  print  intervals  passed? 
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<£9 


Subroutine  INTEG  -  (concluded) 


Does  Error  size 
require  restarting? 


Correct  next 
values  of  0',  0,  0" 


Calculate  error 
estimate:  difference 
between  0"  pred.,  0"  corr. 


Would  doubling  the  step  size 
cause  end  of  print  interval  to  fall 
inside  a  computing  interval? 


Set  Indicator  to  "Double  Step  Size", 
double  computing  step 


Subroutine  DIFFC 


^rn'esm  om&sw  r^nmua  *.««*  ^ 


32 


33 

Subroutine  SCALE  -  (concluded) 


66 


Convert  Argument- 
Card-Field-  Number 
to  fixed  point 


67 


Move  Scale  Table  to 
Punching  region,  convert 
argument  values  to  degrees 
and  write  on  output  tape 


Is  Card  Deck  ready  for  punching? 


3 


@ 


Convert  Deck  ID  to  fixed  point. 
Scale  all  Card  fields  for  punching 
Generate  numbers  for  Calibration- 
Check-Cards, 

spoil  scale  tables  to  prevent 
accidental  re-use 
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APPENDIX  B 

PROGRAM  DESCRIPTION  -  ATMOSPHERE  SUBROUTINE 
NORTH  AMERICAN  AVIATION,  INC.  COLUMBUS 
ENGINEERING  PROGRAM  DESCRIPTION 


1.  Identification 


a.  Atmosphere  Subroutine,  1F113 

b.  Sarah  Crooke,  August  1959 

c.  Aero  Advanced  Design  Group  -  2l4 

d.  Fortran  Source  Deck  is  up  to  date. 

2.  Purpose 

Calculates  standard  atmosphere  properties  at  a  specified  altitude. 
3-  Restrictions 


Fortran  2  subroutine  subprogram. 

4.  Method 


Table  look-up  of  properties  at  base  geopotentials  followed  by  applica¬ 
tion  of  the  ideal  gas  law. 

5.  Use 

Subroutine  is  reached  by  Fortran  Statement 

CALL  A0*$S  (H,  TBL,  R) 

where  H  is  the  input  altitude 

TBL  is  the  table  of  base  geopotentials  and  properties 
R  is  the  results  array. 

6.  Coding  Information 

a.  This  routine  occupies  309  Locations. 

b.  Inferior  Subroutines  are  Library  routines  DUMP,  EXP,  SQRT,  and  the 
Fortran  system  routine  EXP(3* 

7.  Data  Arrangement 

a.  The  array  TBL  must  contain  the  following  data: 


Location 

Description 

Symbol 

1 

Sea  Level  Temperature,  °R 

T. 

2 

Sea  Level  Pressure,  lb/ ft2 

P. 

52 


Location 

Description 

Symbol 

3 

Sea  Level  Gas  Constant,  ft^/sec^°R 

*o 

4 

Sea  Level  Density,  slugs/ft^ 

Pc 

5 

Number  of  Base  Geopotentlals  in  Table 

K 

6 

» 

- 

* 

Base  Geopotentials,  ft 

• 

• 

5  +  a 

- 

6  +  a 

• 

< 

p 

2  2o 

Base  Gas  Constants,  ft  / sec  R 

• 

5  +  2  a 

- 

6  +  2  a 

• 

* 

0 

Base  Temperatures,  R 

Te 

• 

5  +  3a 

« 

6  +  3  a 

« 

* 

Base  Specific  Heat  Ratios,  dimension¬ 

• 

• 

less 

5  +  4a 

% 

6  +  4a 

v 

* 

&  * 

A i*IO  ‘is  the  ith  Base  Pressure 

7  +  4a 

h 

Ratio,  dimensionless 

6ft 

4  +  6a  A* 

5  +  6*  B  J 

b.  On  return  the  array  P  vill  contain  the  following  data: 

Symbol 

4i 
<r 
a 

A 

6 

5 


Location 

1 

2 

3 

4 

5 

6 


Description 

Geopotential  Altitude,  ft 

Density  Ratio,  dimensionless 

Speed  of  Sound,  ft/sec 

Speed  of  Sound,  knots 

Pressure  Ratio,  dimensionless 

Local  Acceleration  of  Gravity, 
ft/sec^ 


*=sa=53ffisigff 


Location  Description  Symbol 

7  Temperature,  °R  T 

8  Gee  Constant,  ft^/ sec^°R  fi. 

9  Incompressible  Dynamic  Pressure 

at  Mach  1,  Ib/ft^  G 

10  Ratio  of  Specific  Heats, 

dimensionless  y 

Formulas 


In  the  formulas  below  H  is  the  input  altitude.  The  remaining  symbols 
are  those  listed  in  7  above.  The  subscript  &  refers  to  the  value  of  the 
quantity  at  the  ith  base  geopotential  altitude,  where 


L  *  I  if  A  <  lx  , 

i  *  j  if  ij  -  ^  fij+i  j  j  *  2  j  •••,  *i  -  /  , 
Ls  n-/  if  An  s 


The  prefix  A  denotes  the  increment  of  the  quantity  from  the  ith  to  the 
( i  ♦/  )th  base  geopotential  altitude. 


A  *  10  811  POOH 

H  *  20  891  000 

T  *  Tt+  (4-  A$)  A  T 

IT 


A  s  A«  + 


6  5 


A  g  T 
A  7* . 


IT 

32. i 7*-  0M  A*  at> 


i/  T*  a  A  -  R^aT  +  0 

r  -3Z./7V  0^9  A%)/aA 

e  > 

if  T*  A  A  -  A*aT*D 


y  *  +  ( l  -  A* )  a  / 


3 


=  51.  <74-0  4-  9 


2  0  8<?/  ooo 


/-/+  20  891  000\ 


a  - 


A  - 


1.6878099 


6  *  I  (a  cr)dl 
Z 
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ii 
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f 


55 


o  in 

o 

m 

o 

■n  o 

n 

o 

n 

o 

in 

o 

in  o  m  o 

m 

o  n 

o 

m  o 

in 

O  in 

o  in 

O 

in 

o 

m  o 

m 

o 

m 

o 

in  o 

n 

o 

o  o 

rH 

r-4 

(M 

c\i  co 

co 

4 

4 

m 

m 

vO 

f-  r-  qo 

ao 

O'  O' 

o 

o 

rH 

rH 

CM  CM 

co  cn 

4 

4 

in 

in  vO 

vO 

r- 

Is- 

CO 

oo  co 

On 

o 

H  H 

H 

H 

H 

H 

rH 

rH 

H 

H 

rH 

H 

H 

rH 

H  H  H 

rH 

rH  rH 

CM 

rvj 

CM 

OsJ 

CM  CM 

CM  CM 

CM 

CM 

CM 

CM  CM 

CM 

CM 

CM 

CM 

(Nl  (M 

CM 

co 

o  o 

o  o 

o 

o 

o  o  o 

O 

o 

o 

o 

o 

o  o  o  o 

c>  o 

O 

o 

o  o 

o  o 

o  o  o 

o 

O 

o  o 

O 

o 

O 

o 

o  o 

o 

o 

4  4 

4 

4 

4 

4 

4- 

4 

4 

4 

4- 

4 

4 

4 

4  4  4 

4 

4  4  4 

4 

4 

4 

4  4 

4  4 

4 

4 

4 

4  4 

4 

4 

4 

4 

4  4 

4 

4 

o  o 

o  o 

o 

o  o 

o 

o 

O 

o  o 

O 

O 

o  o  o  o 

o  o 

o  o 

o 

O 

O  O 

O  O  O 

o 

O 

o  o 

O 

O 

O 

O 

o  o 

O 

O 

u.  u_ 

u. 

u.  u. 

u_ 

u. 

U- 

U. 

u. 

u. 

•u 

u_ 

u.  u.  a.  u.  u. 

U-  U_ 

li-  u. 

u. 

Ll. 

li.  U.  U.U. 

U. 

u. 

u. 

U.  Ul 

li. 

li. 

u_ 

li. 

u.  u. 

u. 

li. 

H  H 

H 

rH 

H 

H 

H 

H 

H 

rH 

H 

rH 

H 

rH  rH  rH 

rH 

rH  rH 

rH 

rH 

rH 

rH 

rH  rH 

rH  rH 

rH 

rH 

rH 

rH  H 

rH 

rH 

rH 

rH 

H  rH 

rH 

rH 

(M 

cm  — 

♦  oc 

*  * 


. — . 

4 

a: 

X 

N 

(30 

+ 

UJ 

ON 

— 

• 

-1 

00 

4 

a: 

o 

CD 

w 

O 

»}c 

o 

< 

o 

oc 

4 

o 

K- 

rH 

♦ 

f— 

Is- 

rH 

w 

«T* 

rH 

O' 

\ 

r** 

QC 

— * 

cn 

r- 

• 

oc 

CO 

— » 

CM 

— »  r 

• 

+ 

■««r 

CM 

■«— 

o 

CM 

— ■ 

07  ~ 

o 

”) 

a 

cn 

CM 

— ■ 

oc 

UJ  o 

o 

n 

— • 

>— 

«— 

- — . 

— ■ 

UJ 

* 

or  rH 

o 

rH 

<n 

ID 

\ 

♦ 

CO 

s 

3 

— > 

» 

rH 

r 

-J 

— 

* 

' — 

• 

CD 

4 

LU  l/) 

0* 

in 

co 

CO 

o: 

CD 

O' 

UJ 

o 

< 

• — 

_J  u 

QO 

rH 

— «■ 

*— 

\ 

< 

*— 

— i • 

-1 

o 

h- 

LiJ 

CD  oc 

o 

r 

o: 

a; 

h- 

oc 

CD 

CD 

o 

_l 

<  • 

CM 

o 

\ 

CD 

♦ 

X 

♦ 

f- 

< 

rH 

— « 

CD 

►— 

* 

CM 

*— 

■r-» 

2: 

* 

o 

O' 

* 

h- 

O' 

in 

< 

•  o 

r 

— » 

CD 

CD 

< 

• 

h~ 

4 

★ 

ON 

CO 

— 

On  i~ 

X  o 

— * 

CD 

a: 

f- 

o 

O 

a 

O 

CO 

o 

QC 

On  * 

rH 

s 

rH 

X 

l 

1 

i 

rH 

• — . 

■ — 

♦ 

o 

4 

• — 

rH 

OC 

CM 

* 

o 

O)  *- 

• 

— 

| 

* 

On 

O' 

cd  cn 

r- 

QC 

•»r 

* 

■W 

00 

<n 

O  uj 

o 

QC 

— » 

in 

m 

m 

— 

■WH 

oc 

•* 

rH 

UJ 

X 

* 

o 

r- 

X  -1 

o 

1 

in 

+ 

A. 

4 

CM 

a: 

a: 

0 

in 

• 

v. 

o 

-J 

a 

O' 

CO  oc 

tD 

o 

O 

— * 

+ 

— ** 

3 

— 

3 

+ 

"S 

+ 

♦ 

* 

X  CM 

CM 

a 

CD 

o 

4 

W • 

NO 

* 

o 

<  < 

rH 

rH 

rH 

m 

4 

4 

w 

4 

~> 

~> 

CD 

X 

X 

o 

* 

cn 

r- 

CD 

< 

Q  O 

a: 

• 

* 

h- 

O' 

•» 

— 

+ 

+ 

•— 

+ 

UJ 

■f 

UJ 

UJ 

W 

X 

o 

a 

a  o 

1 

a 

t- 

+ 

4 

—» 

rH 

CO 

— 

o 

ID 

CO 

in 

u_ 

VC 

UJ 

~> 

-J 

-J 

UJ 

UJ 

1 

I- 

oc 

o  n 

oc 

* 

* 

CD 

r- 

u. 

N. 

— 

o 

r 

Z  Z 

O 

in 

r 

a 

_j 

CD 

> 

CD 

-J 

CD 

_J 

o 

a 

♦ 

— * 

Ll 

* 

X 

rH 

L- 

— * 

cc 

rH 

«— 

o 

*-  o 

CM 

w 

in 

V 

♦ 

ID 

ID 

ID 

CD 

ID  < 

UJ 

< 

CD 

< 

CD 

rH 

■f 

♦ 

CD 

a 

CD 

CM 

in 

< 

• 

a: 

™  * 

• 

a 

* 

*-  »-« 

* 

UJ 

3 

CM 

-J 

3 

-J 

< 

-1 

f- 

_J 

1- 

< 

J— 

v: 

< 

1 — 

CD 

£0  1-0 

X 

a: 

w 

e> 

CM 

o 

m 

H 

II 

o 

3  in 

X 

-J 

rH 

Q 

tl 

CD 

Z 

CD 

h- 

CO 

•W 

CD 

- — 

K 

•*-» 

■f 

1— 

u: 

I- 

OC 

II 

rH 

UJ 

m 

■*— 

VC 

oc 

ti 

CO 

m 

Of 

• 

II 

r«. 

Z 

•> 

o  z 

II 

CD 

1 

M 

X  < 

»— *■ 

< 

II 

v  < 

II 

H* 

< 

II 

*— • 

II 

II 

»— t 

II 

ii 

II 

II 

— 

*— 

ii 

cn 

ii 

n 

ii 

II 

ii 

II 

ii 

M 

»-• 

QC  O 

OC  ID 

r» 

< 

—* ■ 

o  ♦- 

I- 

1- 

— • 

•f 

K— 

r  + 

1- 

X 

CD 

X 

+ 

O 

— * 

— * 

o  a: 

O 

*-• 

**• 

•r-» 

o 

•— 

r-* 

— * 

o 

3 

W' 

CDX 

rH 

H- 

w 

rH 

CM 

— * 

Z 

M 

M 

co  ♦-» 

H 

o 

M 

a 

n 

Z 

o 

•> 

a 

On 

CO 

rH 

W 

CM 

i- 

CM 

in 

CM 

rH 

vD 

4  O' 

4-  io 

1— 

Q 

3  -« 

w 

H 

u. 

— 

O  li. 

o 

H 

CO 

ii 

CO  oc 

M 

D 

»- 

M 

< 

o 

N 

CD 

*— 

w 

li. 

•— 

o 

w 

*— 

W 

— * 

w 

W 

W 

O  uj 

UJ 

z 

to  o 

oc 

*-« 

oc 

Q 

V 

vr  x 

oc  ->  a 

OIK 

Q  “)  O 

o 

3 

a 

OC  OC 

OC 

a: 

OC 

e) 

a: 

a 

oc 

QC  OC 

a 

oc  oc 

oa: 

cr 

ID 

O 

o 

m 

m 

o  m 

o 

m  *h 

CM 

rH 

CM 

cn  <o 

4 

*— « 

H  r— < 

♦— 

U  4-  o 

< 

O  *-« 

u 

— J  o  o 

o 

o  o 

-J 

-*  o 

_J 

ZH  ^ 

< 

u_ 

►— 

u 

o 

in  o 

o 

Z  ^  v* 

Z 

u. 

< 

UJ 

u> 

o: 

UJ 

CO 

51 

3 

z 

(\J  f\J 

< 

ucn  h 

-1 

O  h  rn 

— J  o  o 

51 

o  o 

o: 

o 

O  CY 

Ll_ 

Z  CO 

u. 

*—» 

< 

2 

C£ 

o  m 

UJ 

z  oj  co 

♦— 

Ul 

z 

UJ 

o 

z 

•— * 

Q 

z 

o 

CL 

co  cvj 

CO 

o  h  r- 

UJ 

O  »H  f\J 

QC 

-J  o  o 

QC 

o  o 

O 

U 

r-  r-i 

Z  CO 

X 

iu 

f- 

M 

}C 

o  o 

i/> 

Z  ^  CO 

QC 

u. 

UJ 

UJ 

'O  <t 

uom 

O  H  f\| 

-J  o  o 
o  o 

vO  on 

Z  eg 


• —  o 

j— 

u  o 

a 

u 

o  o 

z 

o 

o 

< 

CO 

z 

u 

o 

u 

UJ 

►— « 

UJ 

Q 

n/> 

z 

UJ 

Q 

CL 

z 

V 

b 

Q 

O 

z 

o 

z 

cc 

< 

H 

UJ 

*“  o 

CL 

fr- 

u  o 

O. 

u 

O  O 

o 

< 

U) 

o 

rH 

UJ 

u 

-J 

u 

UJ 

CQ 

UJ 

Q 

< 

M 

a: 

o 

Q. 

< 

>< 

UJ 

TlGNS  FOR  V 

eg 

< 

►-  o 

u 

1-  s0 

u  o 

o 

u 

o  o 

-1 

O  O 

o 

UJ 

o 

rg 

o 

00 

U 

< 

u  o 

UJ 

a 

UJ  CO 

o 

o 

fr¬ 

iz) 

o 

co 

QC 

Is-  m 
1-4  0 
UOO 

o  o  o 
o  o 


nO  cn 
K  m  ■— c 
O  O'  -O 

o  o  o 
o  o 


m  i-h 

g>  h- 

1  O'  ON 

u  oo  g> 

1  rg  CM 

2: 

UJ  cm  rg 

Q 

X 

< 

O 

QC 

fr— 

O 

H 

a 

O 

o 

(X 

CL 

KJ  CO 

>- 

Ul 

QC 

U 

< 

QC 

CC 

Z> 

CO 

o 

*— • 

0^0 

h-  — c 

-J 

in  <t  4- 

1“  CO  rH 

O  o  o  o 
o  o  o 

vO  CM  00 
U  CMO  CO 

LU  cm  cm  CM 
O 


CD  ->  CD 

2:  >~ 

< 

13 


UOO 

o  o  o 
o  o 

r—  in 
u  co  ■D 
UJ  CM  CM 
Q 


Hinn 
b-  in  o  O 
u  o  o  o 
o  o  o  o 
o  o  o 

l"~  C*>  O' 
U  O'  On  CO 
CM  CM  CM 
Q 


Olo 
x  o  a: 
o 
o 


<M  nO  CM 

o-  in  «o  -0 
u  o  o  o 
o  o  o  o 
o  o  o 

®oo 

o  O'  o>  o' 

UJ  CM  CM  CM 
Q 


*->  X  o 

O  CL 
CC  CO 
Q  Q. 


NO  CO  CM 

t—  cm  cn  cn 
U  ^  -O'  -0 
O  o  o  o 
o  o  o 

no  m  cm 
U  00  CO  CO 
UJ  CM  CM 

Q 


cm  in  o 
»-  in  CM  CM 
U  CM  O  O 
O  O  o  o 
o  o  o 

ON  IM 
U  f-  N-  f- 
UJ  f-i  CM  CM 

o 


■ssesaaas; 


57 


APPENDIX  C 

PROGRAM  DESCRIPTION  -  FIXED  DECIMAL  CARD  READ 
ROUTINE  WITH  ERROR  DETECTION 
NORTH  AMERICAN  AVIATION,  INC.  COLUMBUS 
ENGINEERING  PROGRAM  DESCRIPTION 


1.  Identification 

a.  Fixed  Decimal  Card  Read  Routine  with  Error  Detection  -  Deck  NCF-FXRCH 

b.  0.  C.  Juelich,  March  i960 

c.  Aero  Advanced  Design  Group  -  2l4 

d.  Fortran  Source  Deck  is  up  to  date. 

2.  Purpose 

a.  Load  Input  Data  for  Fortran  II  programs. 

b.  Permit  resumption  of  loading  after  cards  with  zero  or  fractional 
origin  are  found. 

c.  List  such  error  cards  in  the  print-out. 

3.  Restrictions 

a.  Reads  input  tape  5>  writes  output  tape  6  in  case  of  error. 

b.  Data  cards  are  written  on  Form  114-C-17  (Fortran  Fixed  10  Digit 
Decimal  Data) . 

c.  Uses  Fortran  II  Error  Procedure. 

4.  Method 

a.  Card  images  are  read  from  the  input  tape.  The  first  word  specifies 
the  location  to  which  the  data  is  to  be  transmitted,  and  whether  the 
data  set  is  complete. 

b.  Filled- in  card  fields  are  transmitted  to  the  calling  program,  blank 
fields  are  not  transmitted. 

c.  Error  Card  images  are  written  on  the  output  tape.  If  error  cards  are 
found  the  error  procedure  is  invoked  at  the  end  of  the  data  set. 

5.  Use 

a.  Subroutine  is  reached  by  the  Fortran  Statement: 

CALI.  FXRCH  (D)  . 

b.  A  card  image  is  read.  Ihe  first  word  on  the  card  is  the  subscript  h  , 
written  with  decimal  point.  The  ( i  +  Z  )th  word  on  the  card;  1=0, 

1,  2,  3>  4;  is  transmitted  to  location  DHn/  +  i )  as  a  floating 
point  number  unless  it  is  blank  or  -0.  If  the  (  i  +1  )th  word  is 
blank  or  -0  location  BHal*  i )  is  unmodified.  If  n  is  positive  this 
procedure  is  repeated.  If  a  is  negative  control  returns  to  the 
calling  program.  (Subscripting  as  in  BQUTVAIENCE  statements.) 
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6.  Coding  Information 

a.  This  routine  requires  17^  locations. 

b.  Inferior  subroutines  are  the  Library  routine  EHR0R  and  Fortran  System 
Subroutines. 


i 

3 


! 


| 
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APPENDIX  D 

PROGRAM  DESCRIPTION  -  PUNCH  SCALING  HJNCTION 
NORTH  AMERICAN  AVIATION,  INC.  COLUMHJS 
ENGINEERING  PROGRAM  DESCRIPTION 


1.  Identification 

a.  Punch  Scaling  Function,  Deck  -  NCF  NPSCL. 

b.  J.  B.  Burnett,  May  i960. 

c.  Digital  Computing  Group  -  331* 

d.  Source  Deck  is  up  to  date. 

2.  Purpose  ' 

Converts  floating  point  data  to  fixed  point  data  which  is  scaled 
suitably  for  automatic  plotting  equipment. 

i 

3.  Restrictions 

Operates  as  a  Fortran  II  Function  Subprogram. 
k.  Method 

Linear  Scaling .  4 

5.  Use 

i , 

a.  Reached  through  a  Fortran  Statement  such  as 

IFNCH  =  NPSCL  (A,  S) 

where  A  is  the  umber  to  be  scaled 

S  is  an  array  of  dimension  4  containing  the  scale  table. 

S(2)  holds  the  plotter  counts  corresponding  to  the  left 
or  bottom  edge  of  the  graph 

S(4)  holds  the  plotter  counts  corresponding  to  the  right 

or  top  edge  of  the  graph  j; 

S(l)  holds  the  value  of  A  corresponding  to  S(2)  j 

S(3)  holds  the  value  of  A  corresponding  to  S(k)  | 

b.  S(2)  and  S(4)  must  be  between  0.0  and  959*0,  j 

S(3)  must  exceed  S(l).  If  the  array  S  does  not  meet  these  specifics-  j 

tions  IFNCH  will  be  set  to  -0. 

c.  If  A  is  less  than  S(l),  IPNCH  will  be  set  to  0.  f 

If  A  exceeds  S(3)  IPNCH  will  be  set  to  999*  j 

(The  user  must  make  provision  to  confine  the  plotting  head  to  the 

plotting  bed.)  S 

I 
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d.  It  is  expected  that  the  plotting  machine  be  set  to  20  counts  per 
centimeter,  with  (0,  0)  counts  at  the  lower  left  hand  corner  of 
the  grid. 

6.  Coding  Information 

a.  This  routine  occupies  119  locations. 

b.  There  are  no  inferior  subroutines. 


INTERNAL  FORMULA  NUMBERS  AND  OCTAL  LOCAT IONS 
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APPENDIX  E 


PROGRAM  DESCRIPTION  -  OUTPUT  TAPE  WRITE  ROUTINE 
NORTH  AMERICAN  AVIATION,  INC.  COLUMBUS 
ENGINEERING  PROGRAM  DESCRIPTION 


Ident if  icat ion 

a.  Output  Tape  Write  Routine,  Deck  -  NCF  WRITE. 
t>.  0.  C.  Juelieh,  June  i960. 

c.  Research  Group  -  330. 

d.  Fortran  Source  Deck  is  up  to  date. 

Purpose 

Writes  up  to  ten  floating  point  numbers  per  line  retaining  seven 
digit  accuracy  for  most  numbers. 

Restrictions 


a.  Writes  output  tape  6. 

b.  Intended  primarily  for  numbers  between  10" 3  and  10  7  in  magnitude. 
Seven  digit  accuracy  is  retained  for  numbers  between  1  and  107. 

c.  Source  Deck  includes  70k  SAP  instructions. 

Method 


A  Format  Statement  is  generated,  using  F  conversion  for  numbers  less  than 
10 7  and  E  conversion  for  numbers  107  or  larger  in  magnitude. 

Use 


Reached  through  a  Fortran  Statement  such  as 

CALL  WRITE  (6HiIDENT,  A,  N) 

where  i  is  the  pre-print  space  control  character 

IDENT  is  the  indicative  word  for  the  first  line  of  print. 

(The  second  line,  if  any,  will  be  identified  by  an  11, 
the  third  by  21,  etc.) 

A  is  the  first  word  to  be  printed. 

N  is  the  number  of  words  to  be  printed,  N  >  0 . 

Coding  Information 

a.  This  routine  requires  171  locations. 

b.  All  inferior  subroutines  are  part  of  the  Fortran  System. 
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APPENDIX  F 

SUMMARY  OF  LIBRARY  ROUTINES 


SQRT  Evaluates  the  square  root  of  a  non- negative  floating-point 

number.  Invokes  ERR0R  if  argument  is  less  than  zero. 

C0S  Evaluates  the  cosine  of  an  angle  expressed  in  radians  in  float¬ 

ing-point  form. 

Evaluates  the  sine  of  an  angle  expressed  in  radians  in  float¬ 
ing-point  form. 

Evaluates  the  exponential  function  (  e*  )  of  a  floating-point 
number  less  than  88.  Invokes  ERR0R  is  the  argument  exceeds  88. 

a.  Produces  a  "back-trace"  through  the  program,  listing  the 
point  and  subroutine  at  which  the  error  was  detected,  the 
point  at  which  the  subroutine  was  called,  up  to  the  main 
program . 

b.  Turns  on  the  Divide  Check  indicator  and  tests  the  first 
instruction  of  the  main  program.  If  this  is  "IF  DIVIDE 
CHECK"  control  goes  to  this  statement,  otherwise  DUMP  is 
invoked. 

CHAIN  Loads  a  program  identified  by  Record  and  Tape  Number  from  Tape 

to  Core  and  gives  it  control. 

DUMP  Gives  a  full  Memory  Print-out  and  invokes  EXIT. 

Returns  control  to  the  Monitor  to  end  the  job.  (This  function 
is  also  performed  by  the  inprt  tape  reading  system  routine  if 
an  "End  of  File"  is  encountered  on  the  input  tape.) 
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APPENDIX  G 


SAMPLE  INPUT  AND  OUTPUT 


FORTRAN  FIXED  10  DIGIT  DECIM 


FORTRAN  FIXED  10  DIGIT  DECIMAL  DATA 


FORTRAN  FIXED  10  DIGIT  DECIMAL  DATA 


FORTRAN  FIXED  10  DIGIT  DECIMAL  DATA 
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